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Abstract 

We present the Monte Carlo with Absorbing Markov Chains (MCAMC) method for extremely 
long kinetic Monte Carlo simulations. The MCAMC algorithm does not modify the system dynam- 
ics. It is extremely useful for models with discrete state spaces when low-temperature simulations 
are desired. To illustrate the strengths and limitations of this algorithm we introduce a simple 
model involving random walkers on an energy landscape. This simple model has some of the 
characteristics of protein folding and could also be experimentally realizable in domain motion in 
nanoscale magnets. We find that even the simplest MCAMC algorithm can speed up calculations 
by many orders of magnitude. More complicated MCAMC simulations can gain further increases 
in speed by orders of magnitude. 
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I. INTRODUCTION 



There are many excellent algorithms to decrease the computer time required to perform 
Monte Carlo simulations for the statics of a model system. For example, see the articles 
in this volume by D.P. Landau and by B. Diinweg. These algorithms and most other 
acceleration algorithms change the underlying dynamics of the system, which is permissible 
if and only if just the statics of the model is of physical interest. 

However, sometimes the dynamics of the simulation is physically relevant. One example 
is metastability in the Ising model 0, || . By coupling a lattice of quantum spin | particles 
to a quantum heat bath it is possible to obtain the time-dependent density matrix of the 
quantum system of particles plus the bath. If the mixing of the quantum bath is much 
faster than the relaxation of the spins (and using a few other assumptions), the quantum 
bath can be integrated over to obtain a time- dependent Master equation || [|, |[ for a 
classical spin | Ising model on the same lattice. The transition rates in the Master equation 
are related to physical constants such as wave velocities as well as to expectation values in 
the original quantum system. This calculation was performed with a particular fermionic 
bath by Martin in 1977 0, and the Glauber dynamic H was obtained. Recently, Park et 
al. 0, |U performed the same calculation with a (^-dimensional bosonic bath and obtained 
somewhat different transition rates in the Master equation. This bosonic dynamic is relevant 
for molecular magnets, and also leads to some novel features in metastable decay M. The 
dynamic obtained in the above fashion is physical, and hence cannot be modified if the 
time-dependence of the model is to be compared with experiments. Note that this approach 
relates the dynamic Monte Carlo simulation time, measured in Monte Carlo Steps per Spin 
(MCSS), to the physical time (in seconds). 

For the Ising simulation, the attempt frequency in the kinetic Monte Carlo is related to 
the inverse phonon frequency, about 10~ 13 s. To study metastable decay the time scales are 
typically on the order of human times (seconds to years), or for relevance to paleomagnetism 
the time scale is many millions of years. Consequently, algorithms which do not change the 
dynamic but are faster-than-real-time must be used to directly compare with experiments. 
This paper details one such algorithm, the Monte Carlo with Absorbing Markov Chains 
(MCAMC) algorithm H, 110]. A recent review at the introductory level of faster-than-real- 



time kinetic Monte Carlo algorithms, including the MCAMC algorithm, is available [|TT 



II. MODEL 



MCAMC simulations for the Ising model have been presented previously || [10], [11], |12 
They are related to magnetization switching in thin nanoscale highly-anisotropic magnetic 
films. In this paper we introduce a simple model to illustrate the MCAMC method. We will 
see below that this simple model also has some interesting physics, and also can be related 
both to magnetization switching of thin films and to questions related to protein folding. 

Consider a one-dimensional lattice where each site % has been assigned an energy Ei. This 
model is generalizable to higher dimensions, and the MCAMC method will work in higher 
dimensions, but here for simplicity we focus on the linear lattice. In particular, consider the 
20-site lattice shown in Fig. 1, with energies given in Table 1. 

On this lattice we initially randomly place a number of walkers, N w . Introduce the dy- 
namic that at each Monte Carlo step a walker is randomly picked (with uniform probability) , 
and then whether the walker attempts to move left or right is randomly chosen. The chosen 
walker then moves to the adjacent lattice site with a probability 

exp(-^ ±1 /fc B T) 

Pmove ~ expi-Ei/kaT) + exp(- E t±1 / k B T) ' 1 ' 

where the + (— ) sign is for the attempted move to the right (left). We place reflecting walls 

(infinite E) at sites and 21. One attempted move is a Monte Carlo step (mcs), and iV w 
attempts is a Monte Carlo step per walker (MCSW). We will be interested in the average 
time, (r), for all N w walkers to first reach the same lattice site. This 'coagulation' of walkers 
may or may not be at the global energy minima (site 9). 

This model can be viewed as a model for switching in a nanomagnet where a domain 
wall is constrained to lie in a thin film that is anchored at the two ends by macroscopic 
magnets. The energy at each site would then correspond to the energy the domain wall 
would have when it is in a particular coarse-grained location in the thin film. The lowest 
energy points would be where pinning of the magnetic domain wall is the strongest and the 
highest energies would be the saddle points a domain wall would need to thermally overcome 
to traverse from one pinning site to another. The physical question is when N w independent 
thin-film devices which start in a random state (demagnetized) would all first be in the same 
location. 

The same model can be viewed as a very simple example for protein folding by considering 



each point to be the free-energy at some abstract point in phase space |L3] . Then the lowest 



energy corresponds to the native configuration. The physical question would then be the 
average speed at which the protein folds, i.e. when N w proteins would be in the same 
configuration. 



III. KINETIC MONTE CARLO 

A number, M, of different coagulations of walkers will be performed to obtain the average 
lifetime (r), i.e. the average time until coagulation occurs. For each coagulation the iV w 
walkers are first placed randomly, with uniform probability on each site. Then the kinetic 
Monte Carlo simulation is performed, calculating the number of time steps n for coagulation 
i until all N w walkers are simultaneously at the same lattice site. This lattice site may be 
any of the 20 sites of the lattice. The average lifetime is then given by 



The kinetic Monte Carlo algorithm is simple. Each move requires 3 uniformly distributed 
random numbers, Tj, with < < 1. The first random number is used to select uniformly 
one of the N w walkers. This is accomplished by choosing walker j, given by j = 1 + [^i-A^J 



otherwise it attempts a move to the right. Finally if r% < p movc with p move from Eq. ([[]) the 
walker moves to the chosen new lattice position. Whether or not a move is made, the time 
is advanced by one unit, = + 1. 

The results for (r) using this algorithm are shown in Fig. § and || The data will be 
discussed in the results section, but here and in the next two sections we concentrate on the 
algorithmic aspects. The CPU (central processing unit) time required (on a single processor 
of a Cray SV1 vector computer) for this algorithm is shown in Fig. f| with the label mcl. 
Note that at both very low and very high temperatures this algorithm requires substantial 
amounts of computer time. In particular, the average CPU time required for this algorithm 
is proportional to the value of (r) . Thus at low temperatures where (r) grows exponentially 
fast in T _1 the time required for the simulation also grows exponentially quickly. 




(2) 



where |_^J is the integer part of x. If r<i < -, the chosen walker attempts a move to the left, 



IV. n-FOLD WAY = s=l MCAMC 



To decrease the simulation time required, without altering the dynamic, an event-driven 
simulation can be performed. This is also called an ra-fold way simulation |TJJ . The original 



paper |TJ| used continuous time, but the same algorithm can be cast into the discrete-time 
version (where Tj = Tj + 1 at each time step) (T5| used here. The n-fold way algorithm is 
a s=l MCAMC algorithm, because the system has one transient state (s=l), which is the 
current state of the system. In our case it has 2N W absorbing states (states the system 
can jump to from the current state), two (one of which may have zero probability) for each 
random walker. 

The n-fold way algorithm also requires 3 random numbers, r« at each step. First form a 
vector Pj ump of length N w which contains the N w probabilities Pj nmp (k) that a walker jumps to 
either the left or right given that it was picked in the uniform picking part of the algorithm. 
Then increment the time (in mcs) by 



Ti = n + 



ln(rx) 



In 1 - 



+ 1 (3) 



with p sum = J2k=i Pjump(k). Next the walker, j, that actually jumped is calculated by finding 
the value of j that satisfies 

j'-i i 

^Pjump(&) < r 2Psum < ^Pjumpik) (4) 
k=l k=l 

where the first sum is taken as zero if j = 1. Finally if the probability of moving left, p move -, 
from Eq. (|l]) satisfies the relation 

Pmove - < r 3 Pj ump (j) (5) 

the j th walker is moved to the left, otherwise it is moved to the right. This algorithm is 
repeated until all walkers coagulate at one point. 

The results for this discrete-time n-fold way algorithm is identical (within statistical 
errors for these M = 10 3 coagulations) to those of the standard kinetic Monte Carlo. The 
way the dynamic has been implemented on the computer is just different. Results from 
this algorithm are labeled nfl in the figures. The n-fold way algorithm requires additional 
calculations at each step, but the time increment that is added to Ti may be larger than unity 



at each step. This can drastically decrease the simulation time required, particularly at low 
temperatures, as seen in Fig. f|. This is because the n-fold way algorithm is an event-driven 
algorithm. In other words the time is incremented only when an event happens, namely 
only when a walker jumps from its current site. The number of time steps before one walker 
jumps can be very large, particularly at low temperatures. As seen in Fig. [| the n-fold way 
algorithm requires about an order of magnitude less CPU time at low temperatures than 
the previous implementation of the algorithm. A factor of 10 is important in simulations, 
but at low temperatures where (r) is growing exponentially with T _1 , so does the required 
simulation time using the n-fold way method. 



V. MCAMC WITH s=2 



At low temperatures it would be very nice to have a faster algorithm than the n-fold way. 
The reason the n-fold way algorithm scales so poorly (exponentially in T _1 ), is because 
when a walker is in a flat energy minima (at sites 5 and 6 or at sites 16 and 17) the average 
time before the walker moves is equal to a value which is independent of temperature. In 
particular, given that the walker in one of these minima is picked it has a probability of ■? 
of moving to the adjacent equal-energy site. Thus at low temperatures the walker in such a 
minima will rattle back and forth many times before it jumps. 

To user higher s MCAMC, more states are included in the transient subspace [fTE ], One 
way of doing this for a model closely related to the current model is given in [IT] . That way 
is easily generalized to cases where the energies in the minima are nearly equal, or to the 
case of including larger numbers of transient states in the calculation. Here we introduce a 
simple method that works most easily when the energies in some minima are equal and the 
energies to hope from these minima are equal. 

Let N s i be the number of walkers that are not located at sites 5-6 or 16-17, and the 
number that are in one of these minima. Clearly N w = N sl + N s2 . If all the walkers are 
located in either the 5-6 or the 16-17 minima, the n-fold way algorithm of the previous section 
must be used. Otherwise, consider the current state of the system to be expanded to be the 
state with the N sl walkers fixed at their current site and the N s2 walkers are still located in 
their respective 5-6 or 16-17 minima. Form the vector of length N w with elements either: 1) 
if the walker is not at the 5-6 or 16-17 minima, the (n-fold way) probability Pj um _ p (k) that 



a walker jumps to either the left or right given that it was picked in the uniform picking 
part of the algorithm; 2) if the walker is at the 5-6 or 16-17 minima, the probability that 
the walker exits this minima given that it was picked during the random picking part of the 
algorithm. Then increment the time (in mcs) by Tj = Tj + At; with 



An 



ln(ri) 



In 1 



Ms 1 1 n 



+ 1 (6) 



with p sum = X]fc=LPjum P (^) an d ri a uniformly distributed random number. 

Next use a random number T2 to pick one of the j walkers to move, in the same fashion 
as was done for the n-fold way method, Eq. |j. The only difference here is that when one 
of the N s2 walkers is picked it will exit from the 5-6 or 16-17 minima, not just move to an 
adjoining lattice site. 

Next use random numbers to find the number of times m ; that each of the iV w walkers 
was picked given that the system exited from the current state at time At*. One way of 
doing this is with a tree-like structure. For example, for iV w = 4 use 3 random numbers 
rs, Ti and r$, with the number of times walker 1 was picked equal to mi = r4r3(ATj — 1), 
the number of times walker 2 was picked m% = (1 — r^r^An — 1), for walker 3 one has 
7713 — r s(l — r 3)(Arj — 1), and for walker 4 one has 777.4 = (1 — r s)(l ~~ r 3)(Arj — 1). Make sure 
that rounding effects does not change At,, i.e. ensure that Ar^ = 1 + Y2f=i m (- For each of 
the N S 2 walkers use its own random number r x , and move walker i to its other equal-energy 
minimum site if 

1 1 / l x 



r*<- 2 -- 2 {- 2 ) • « 

otherwise leave it where it was. 

If the j th walker that was picked to move is in the 5-6 or 16-17 minima, move it out 
from the minima to its adjacent higher-energy site. If the j th walker is not in one of these 
minima, then as in the n-fold way, use a random number r y and the probability of moving 
left, p m ove,-, from Eq. ([[]), and move the j th walker left if 

Pmove,— — ^j/Pjump (jj ) i (8) 



otherwise it is moved to the right. This algorithm is repeated until all walkers coagulate at 
one point. 



As seen in the figures, with this algorithm denoted by amcl , the average lifetimes for this 
algorithm are also equal (within statistical errors for M = 10 3 escapes) to the results from 
the other two algorithms. However, as seen in Fig. |], the time required for the simulation 
at low temperatures is approximately independent of temperature! Thus this algorithm is 
many times faster at low temperatures than even the n-fold way algorithm. At T _1 = 50 
this algorithm is about 10 8 times faster than the n-fold way algorithm. 

The reason this algorithm scales so well at low temperatures is that in one algorithmic 
step one of the walkers moves to a site which has a higher energy, and consequently such a 
move will become less probable at low temperatures. The jump probability of Eq. (|7|) can 
be obtained using a Markov chain with the probability of remaining at the current lattice 
site equal to | and the probability of going to the other equal-energy site equal to 7. Note 
that we already know this walker stays in its minima for rrii time steps, which allows us to 
use a Markov chain. 

In general higher s MCAMC algorithms could also be used. For certain energy landscapes 
(for example where a walker will rattle around in a local energy minimum composed of more 
than two sites) such higher s MCAMC algorithms will be needed to obtain an algorithm 
that scales approximately independently of T for low temperatures. 

VI. RESULTS 

The results in Fig. ^|and |]show that there is a minimum in (r) with temperature. At low 
temperatures for large iV w , most likely there is at least one walker that will be in the 16-17 
minimum, and this walker will require a long time to move over the saddle point (at site 12) 
to join the other walkers (which by that time will probably all be at the global minimum, 
site 9). At high temperatures, the number of steps before all walkers will be at the same 
site should be approximately the same as the probability of them all being at the same site 
if they are placed randomly on the lattice, a probability given here by (^) Nw For high 
temperatures this can be seen in Fig. 

Thus there is a temperature where the lifetime (r) is a minimum. Furthermore, the 
approximate width of this minimum decreases with N w , as seen in Fig. This should be 
expected in other models with many degrees of freedom that have complicated energy for 
free-energy landscapes, such as models for protein folding. 



VII. CONCLUSION AND DISCUSSION 



The Monte Carlo with Absorbing Markov Chains (MCAMC) algorithm was introduced 
and applied to a simple model. This model could be realized experimentally (in a course- 
grained fashion) using nanoscale magnetic films with non-constant widths. It can also be 
viewed as a toy model to study the average time in which a protein will fold, or any other 
such model with a complicated energy surface and intrinsic dynamics. The model shows 
that there is a minimum in the average lifetime (r) as a function of temperature. At low 
temperature the lifetime is large because a large energy barrier must be overcome. At high 
temperatures the lifetime to coagulation is large because it is improbable that the walkers 
will remain in a low energy configuration very long. 

At low temperatures, the MCAMC algorithm gives exponentially fast speed-ups compared 
to the traditional kinetic Monte Carlo algorithm. The s=2 MCAMC algorithm in fact scales 
almost independently of temperature at low temperatures, compared to an exponentially 
growing simulation time required for the traditional kinetic Monte Carlo and the n-fold 
way simulations. Such faster-than-real-time algorithms are required for realistic dynamic 
simulations of many model systems. 
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FIG. 1: The energies of the 20 sites are shown, and are listed in Table 1. 
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FIG. 2: The average lifetime per walker, (t)N^ 1 from 10 3 escapes is shown as a function of T . 
Note the logarithmic scale. The solid lines join the points for 6 walkers. The 9 different points are 
for 3 different values of walkers (iV w = 4, 6, 8) and 3 different programs labeled mcl for normal 
Monte Carlo, nfl for ra-fold way, and amcl for s=2 MCAMC. 




FIG. 3: The same as Fig. 2 but plotted as a function of T to show the behavior at high temp 
tures. All symbols and notation is the same as Fig. 2. 
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FIG. 4: The CPU time required in minutes to run 10 3 escapes on a vector Cray SV1 computer is 
shown as a function of T . Note the logarithmic time scale. The solid lines join the points for 6 
walkers. The 9 different points are for 3 different values of walkers (iV w = 4, 6, 8) and 3 different 
programs labeled mcl for normal Monte Carlo, nfl for n-fold way, and amcl for s=2 MCAMC. 



Tables 



TABLE I: The 20 energies used in the simulation, as in Fig. 1. 







Site Number 


Energy 


1 


1.0 


2 


1.0 


3 


0.0 


4 


0.5 


5 


0.0 


6 


0.0 


7 


0.5 


8 


0.0 


9 


-2.0 


10 


0.5 


11 


0.5 


12 


0.75 


13 


0.0 


14 


0.5 


15 


0.0 


lb 


-1.0 


17 


-1.0 


18 


0.0 


19 


0.0 


20 


-1.0 







